# First replication file for Weidmann & Salehyan, Violence and Ethnic Segregation: A Computational Model Applied to Baghdad, ISQ (2013) 57, 52-64
# This file replicates the results presented in Figures 5, 6 and 7, and in Table 1

# Nils B. Weidmann (nils.weidmann@uni-konstanz.de)

library(doBy)
library(apsrtable)

setwd("/Users/nilsw/OldProjects/Baghdad/R")

# read data files
res.all.1 <- read.table("counterfactual1_countPunishedAttacks.txt", header=T)
res.all.2 <- read.table("counterfactual2_countPunishedAttacks.txt", header=T)
res.all.3 <- read.table("counterfactual3.txt", header=T)


# different policing rates (implemented from the beginning)
lower <- function(x) {quantile(x, probs=0.25, names=F)}
upper <- function(x) {quantile(x, probs=0.75, names=F)}
res <- summaryBy(cumNumAttacks + prophomunits ~pCaptureIns, data=res.all.1, FUN=c(mean))
res.std <- summaryBy(cumNumAttacks + prophomunits ~pCaptureIns, data=res.all.1, FUN=c(lower, upper))

# Figure 6a
pdf("counterfactual_violence.pdf", width=7, height=7)
plot(res$pCaptureIns, res$cumNumAttacks.mean, type="l", lwd=3, xlab="Policing success rate", ylab="Total number of attacks", cex.lab=1.5, ylim=c(0, 4200))
dev.off()

# Figure 6b
pdf("counterfactual_segregation.pdf", width=7, height=7)
plot(res$pCaptureIns, res$prophomunits.mean, type="l", lwd=3, xlab="Policing success rate", ylab="Proportion of homogeneous units", cex.lab=1.5, ylim=c(0.3, 0.65))
dev.off()

# different onsets of policing
res.low <- subset(res.all.2, pCaptureIns==0.1)
res.medium <- subset(res.all.2, pCaptureIns==0.3)
res.high <- subset(res.all.2, pCaptureIns==0.5)
reslow.agg <- summaryBy(cumNumAttacks + prophomunits ~startCaptureAt, data=res.low, FUN=c(mean))
resmedium.agg <- summaryBy(cumNumAttacks + prophomunits ~startCaptureAt, data=res.medium, FUN=c(mean))
reshigh.agg <- summaryBy(cumNumAttacks + prophomunits ~startCaptureAt, data=res.high, FUN=c(mean))

# Figure 7a
pdf("cf_onset_violence.pdf", width=7, height=7)
plot(reslow.agg $startCaptureAt, reslow.agg $cumNumAttacks.mean, type="l", lwd=3, xlab="Policing onset (time step)", ylab="Total number of attacks", cex.lab=1.5, ylim=c(0, 4200))
lines(resmedium.agg$startCaptureAt, resmedium.agg$cumNumAttacks.mean, type="l", lty=2, lwd=3)
lines(reshigh.agg$startCaptureAt, reshigh.agg$cumNumAttacks.mean, type="l", lty=3, lwd=3)
legend("bottomright", legend=c("low (0.1)", "medium (0.3)","high (0.5)"), lwd=3, lty=c(1,2,3), cex=1.5, title="Policing")
dev.off()

# Figure 7b
pdf("cf_onset_segregation.pdf", width=7, height=7)
plot(reslow.agg $startCaptureAt, reslow.agg $prophomunits.mean, type="l", lwd=3, xlab="Policing onset (time step)", ylab="Proportion of homogeneous units", cex.lab=1.5, ylim=c(0.3, 0.65))
lines(resmedium.agg$startCaptureAt, resmedium.agg$prophomunits.mean, type="l", lty=2, lwd=3)
lines(reshigh.agg$startCaptureAt, reshigh.agg$prophomunits.mean, type="l", lty=3, lwd=3)
legend("bottomright", legend=c("low (0.1)", "medium (0.3)","high (0.5)"), lwd=3, lty=c(1,2,3), cex=1.5, title="Policing")
dev.off()

# trace simulation over time
lower <- function(x) {quantile(x, probs=0.20, names=F)}
upper <- function(x) {quantile(x, probs=0.80, names=F)}
trace <- summaryBy(numAttacks + prophomunits ~ time, data=res.all.3, FUN=c(mean, lower, upper))

# Figure 5a
pdf("trace_violence.pdf", width=7, height=7)
plot(trace$time, trace$numAttacks.mean, type="l", lwd=3, ylab="Number of attacks", xlab="Time")
dev.off()

# Figure 5b
pdf("trace_segregation.pdf", width=7, height=7)
plot(trace$time, trace$prophomunits.mean, type="l", lwd=3, ylab="Proportion of homogenous units", xlab="Time", ylim=c(0,1))
lines(trace$time, trace$prophomunits.upper, col="grey", lwd=3)
lines(trace$time, trace$prophomunits.lower, col="grey", lwd=3)
dev.off()

# Regression results in Table 1
res.all.3$prophomunits2 <- res.all.3$prophomunits^2
res.all.3$lnumAttacks <- log(res.all.3$numAttacks+1)
runids <- unique(res.all.3$runId)
lm1 <- lm(lnumAttacks ~ prophomunits + prophomunits2 + as.factor(runId), data=res.all.3)
summary(lm1)
